!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MC&MD with FFT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE mc_md(qs, ps, sample_number )

  use prec
  use md_vv
  use emc
  use global_data
  IMPLICIT NONE
  include 'fftw3.f'


    integer, intent(in) :: sample_number 
    real(r8), dimension(:,:), allocatable :: jautoavg, j1autoavg, rms, rms1
    real(r8), dimension(:,:,:), allocatable :: jauto, j1auto
    real(r8), dimension(:,:), allocatable :: jjarray, jjarray1
    real(r8),dimension(2) :: jtot, j0, j1, j1tot, javg, jjint, jjint1
    real(r8),dimension(2) ,intent(in) :: qs, ps
    real(r8),dimension(2)::us, vs
    real(r8),dimension(2) :: uscopy, vscopy
    integer(i12) :: isample, msample         ! sample index
    integer :: isnap, i
    real(r8) :: emccurrent   ! current value of energy as calculated by emc routines

! for FFTW3
  integer :: FFTW_METHOD
  real(r8), dimension(:), allocatable :: in, in1, in2, in3
  complex(r8), dimension(:), allocatable :: out, out1, out2, out3
  integer*8 :: planforward, planbackward
  integer*8 :: planforward1, planbackward1
  integer*8 :: planforward2, planbackward2
  integer*8 :: planforward3, planbackward3

! for output file
  character (len=20), dimension(2,20) :: correlation

correlation(1,1) = "x_correl1.dat" ! file for autocorrelation fn for delta_j of mode x on trial 1
correlation(1,2) = "x_correl2.dat"
correlation(1,3) = "x_correl3.dat"
correlation(1,4) = "x_correl4.dat"
correlation(1,5) = "x_correl5.dat"
correlation(1,6) = "x_correl6.dat"
correlation(1,7) = "x_correl7.dat"
correlation(1,8) = "x_correl8.dat"
correlation(1,9) = "x_correl9.dat"
correlation(1,10) = "x_correl10.dat"
correlation(1,11) = "x_dcorrel1.dat"  ! file for autocorrelation fn for j_dot of mode x on trial 1
correlation(1,12) = "x_dcorrel2.dat"
correlation(1,13) = "x_dcorrel3.dat"
correlation(1,14) = "x_dcorrel4.dat"
correlation(1,15) = "x_dcorrel5.dat"
correlation(1,16) = "x_dcorrel6.dat"
correlation(1,17) = "x_dcorrel7.dat"
correlation(1,18) = "x_dcorrel8.dat"
correlation(1,19) = "x_dcorrel9.dat"
correlation(1,20) = "x_dcorrel10.dat"
correlation(2,1) = "y_correl1.dat"  ! file for autocorrelation fn for delta_j of mode y on trial 1
correlation(2,2) = "y_correl2.dat"
correlation(2,3) = "y_correl3.dat"
correlation(2,4) = "y_correl4.dat"
correlation(2,5) = "y_correl5.dat"
correlation(2,6) = "y_correl6.dat"
correlation(2,7) = "y_correl7.dat"
correlation(2,8) = "y_correl8.dat"
correlation(2,9) = "y_correl9.dat"
correlation(2,10) = "y_correl10.dat"
correlation(2,11) = "y_dcorrel1.dat"  ! file for autocorrelation fn for j_dot of mode y on trial 1
correlation(2,12) = "y_dcorrel2.dat"
correlation(2,13) = "y_dcorrel3.dat"
correlation(2,14) = "y_dcorrel4.dat"
correlation(2,15) = "y_dcorrel5.dat"
correlation(2,16) = "y_dcorrel6.dat"
correlation(2,17) = "y_dcorrel7.dat"
correlation(2,18) = "y_dcorrel8.dat"
correlation(2,19) = "y_dcorrel9.dat"
correlation(2,20) = "y_dcorrel10.dat"

  us = qs
  vs = ps
  allocate( jjarray(1:2,0:nsnaps) )
  allocate(jauto(1:2,1:10,0:nsnaps))
  allocate(in(0:nsnaps))
  allocate(out(0:(nsnaps+1)/2))
  allocate( jjarray1(1:2,0:nsnaps) )
  allocate(j1auto(1:2,1:10,0:nsnaps))
  allocate(in1(0:nsnaps))
  allocate(out1(0:(nsnaps+1)/2))
  allocate(jautoavg(1:2,0:nsnaps))
  allocate(j1autoavg(1:2,0:nsnaps))
  allocate(in2(0:nsnaps))
  allocate(out2(0:(nsnaps+1)/2))
  allocate(in3(0:nsnaps))
  allocate(out3(0:(nsnaps+1)/2))
  allocate(rms(1:2,0:nsnaps))
  allocate(rms1(1:2,0:nsnaps))

  jautoavg = 0.
  j1autoavg = 0.
  jauto = 0.
  j1auto = 0.
  rms = 0.
  rms1 = 0.

do i=1,10
  call setupemc( du, dv, us, vs, m, seed+i , lambda, emccurrent )
  call setupmdvv( lambda, dt )


  ! initialize  FFTW3
  FFTW_method = FFTW_ESTIMATE
  call dfftw_plan_dft_r2c_1d( planforward, nsnaps, in, out, FFTW_method )
  call dfftw_plan_dft_c2r_1d( planbackward, nsnaps, out, in, FFTW_method )
  call dfftw_plan_dft_r2c_1d( planforward1, nsnaps, in1, out1, FFTW_method )
  call dfftw_plan_dft_c2r_1d( planbackward1, nsnaps, out1, in1, FFTW_method )
  call dfftw_plan_dft_r2c_1d( planforward2, nsnaps, in2, out2, FFTW_method )
  call dfftw_plan_dft_c2r_1d( planbackward2, nsnaps, out2, in2, FFTW_method )
  call dfftw_plan_dft_r2c_1d( planforward3, nsnaps, in3, out3, FFTW_method )
  call dfftw_plan_dft_c2r_1d( planbackward3, nsnaps, out3, in3, FFTW_method )
  
  jauto = 0.
  j1auto = 0.

!+++++++ equilibrate first
!
  call manystepsemc( us, vs, m, nmcequil )
!
!
!+++++++ now start sampling
    
  msample = max(sample_number/100,1)

  javg = 0.
  do isample = 1,sample_number

  	call manystepsemc( us, vs, m, nmcsteps )
		javg(1) = javg(1) + (vs(1)**2 + us(1)**2*omega1**2)/(2*omega1)
		javg(2) = javg(2) + (vs(2)**2 + m**2*us(2)**2*omega2**2)/(2*omega2*m)

  enddo

     javg = javg/(sample_number)

  print *, javg

  do isample = 1,sample_number

     if ( mod(isample-1,msample).eq.0 ) print *,' isample = ',isample

     call manystepsemc( us, vs, m, nmcsteps )

        j0(1) = (vs(1)**2 + us(1)**2*omega1**2)/(2*omega1)
        j0(2) = (vs(2)**2 + m**2*us(2)**2*omega2**2)/(2*omega2*m)
	j1(1) = (vs(1)*us(1)*(1-omega1**2+lambda*us(2)**2))/omega1
	j1(2) = (vs(2)*us(2)*(1-m*omega2**2+lambda*us(1)**2))/(omega2*m)


! copy the initial state so we can also do a trajectory with p ->  -p
     uscopy = us
     vscopy = vs
     
     jjarray = 0.
     jjarray1 = 0.
     jjarray(:,0) = jjarray(:,0) + j0(:)
     jjarray1(:,0) = jjarray1(:,0) + j1(:)
      do isnap=1,nsnaps
         call mdvv( us, vs, m, nmdstepsbetweensnaps )

        jtot(1) = (vs(1)**2 + us(1)**2*omega1**2)/(2*omega1)
        jtot(2) = (vs(2)**2 + m**2*us(2)**2*omega2**2)/(2*omega2*m)
        j1tot(1) = (vs(1)*us(1)*(1-omega1**2+lambda*us(2)**2))/omega1
        j1tot(2) = (vs(2)*us(2)*(1-m*omega2**2+lambda*us(1)**2))/(omega2*m)


         jjarray(:,isnap) = jjarray(:,isnap) + jtot(:)
	 jjarray1(:,isnap) = jjarray1(:,isnap) + j1tot(:)
     enddo


!!!!!fft in the west?!!!!!!!!!!!!!!!!!!!
     in(:) = jjarray(1,:) - javg(1)
     call dfftw_execute(planforward)
     out(:) = abs(out(:))**2
     call dfftw_execute(planbackward)
     in1(:) = jjarray1(1,:)
     call dfftw_execute(planforward1)
     out1(:) = abs(out1(:))**2
     call dfftw_execute(planbackward1)

     in2(:) = jjarray(2,:) - javg(2)
     call dfftw_execute(planforward2)
     out2(:) = abs(out2(:))**2
     call dfftw_execute(planbackward2)
     in3(:) = jjarray1(2,:)
     call dfftw_execute(planforward3)
     out3(:) = abs(out3(:))**2
     call dfftw_execute(planbackward3)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     jauto(1,i,:) = jauto(1,i,:) + in(:)/float(nsnaps)**2
     j1auto(1,i,:) = j1auto(1,i,:) + in1(:)/float(nsnaps)**2
     jauto(2,i,:) = jauto(2,i,:) + in2(:)/float(nsnaps)**2
     j1auto(2,i,:) = j1auto(2,i,:) + in3(:)/float(nsnaps)**2
     
     us = uscopy
     vs = -vscopy

     jjarray = 0.
     jjarray1 = 0.

        j0(1) = (vs(1)**2 + us(1)**2*omega1**2)/(2*omega1)
        j0(2) = (vs(2)**2 + m**2*us(2)**2*omega2**2)/(2*omega2*m)
        j1(1) = (vs(1)*us(1)*(1-omega1**2+lambda*us(2)**2))/omega1
        j1(2) = (vs(2)*us(2)*(1-m*omega2**2+lambda*us(1)**2))/(omega2*m)


     jjarray(:,0) = jjarray(:,0) + j0(:)
     jjarray1(:,0) = jjarray1(:,0) + j1(:)

     do isnap=1,nsnaps
        call mdvv( us, vs, m, nmdstepsbetweensnaps )

        jtot(1) = (vs(1)**2 + us(1)**2*omega1**2)/(2*omega1)
        jtot(2) = (vs(2)**2 + m**2*us(2)**2*omega2**2)/(2*omega2*m)
        j1tot(1) = (vs(1)*us(1)*(1-omega1**2+lambda*us(2)**2))/omega1
        j1tot(2) = (vs(2)*us(2)*(1-m*omega2**2+lambda*us(1)**2))/(omega2*m)

	jjarray(:,isnap) = jjarray(:,isnap) + jtot(:)
 	jjarray1(:,isnap) = jjarray1(:,isnap) + j1tot(:)
     enddo

!!!!!fft in the west?!!!!!!!!!!!!!!!!!!!
     in(:) = jjarray(1,:) - javg(1)
     call dfftw_execute(planforward)
     out(:) = abs(out(:))**2
     call dfftw_execute(planbackward)
     in1(:) = jjarray1(1,:)
     call dfftw_execute(planforward1)
     out1(:) = abs(out1(:))**2
     call dfftw_execute(planbackward1)

     in2(:) = jjarray(2,:) - javg(2)
     call dfftw_execute(planforward2)
     out2(:) = abs(out2(:))**2
     call dfftw_execute(planbackward2)
     in3(:) = jjarray1(2,:)
     call dfftw_execute(planforward3)
     out3(:) = abs(out3(:))**2
     call dfftw_execute(planbackward3)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
     jauto(1,i,:) = jauto(1,i,:) + in(:)/float(nsnaps)**2
     j1auto(1,i,:) = j1auto(1,i,:) + in1(:)/float(nsnaps)**2
     jauto(2,i,:) = jauto(2,i,:) + in2(:)/float(nsnaps)**2
     j1auto(2,i,:) = j1auto(2,i,:) + in3(:)/float(nsnaps)**2    

 
  enddo

! do isnap=0, nsnaps
!  print *, in(isnap)
! enddo

 
 jauto(:,i,:) = jauto(:,i,:)/(2*sample_number)
 j1auto(:,i,:) = j1auto(:,i,:)/(2*sample_number)
 jauto(1,i,:) = jauto(1,i,:)/jauto(1,i,0)
 j1auto(1,i,:) = j1auto(1,i,:)/j1auto(1,i,0)
 jauto(2,i,:) = jauto(2,i,:)/jauto(2,i,0)
 j1auto(2,i,:) = j1auto(2,i,:)/j1auto(2,i,0)
 jautoavg(1,:) = jautoavg(1,:)+jauto(1,i,:)
 j1autoavg(1,:) = j1autoavg(1,:)+j1auto(1,i,:)
 jautoavg(2,:) = jautoavg(2,:)+jauto(2,i,:)
 j1autoavg(2,:) = j1autoavg(22,:)+j1auto(2,i,:)

!  do isnap=0, nsnaps
!  print *, j1auto(1,isnap)
! enddo

 ! write a file with correl fn of delta_j for mode x
  open(file=correlation(1,i),status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9211)isnap*snaptime,jauto(1,i,isnap)
  9211 format(1x,2g24.16)
  enddo
  close(10)

! write a file with correl fn of j_dot for mode x
  open(file=correlation(1,10+i),status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9212)isnap*snaptime,j1auto(1,i,isnap)
9212 format(1x,2g24.16)
  enddo
  close(10)

! write a file with correl fn of delta_j for mode y
  open(file=correlation(2,i),status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9213)isnap*snaptime,jauto(2,i,isnap)
  9213 format(1x,2g24.16)
  enddo
  close(10)

! write a file with correl fn of j_dot for mode y
  open(file=correlation(2,10+i),status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9214)isnap*snaptime,j1auto(2,i,isnap)
9214 format(1x,2g24.16)
  enddo
  close(10)


  enddo
 
  jautoavg(:,:) = jautoavg(:,:)/10
  j1autoavg(:,:) =j1autoavg(:,:)/10

do i=1,10
	rms(1,:) = rms(1,:)+ (jauto(1,i,:)-jautoavg(1,:))**2
	rms1(1,:) = rms1(1,:) +(j1auto(1,i,:)-j1autoavg(1,:))**2
	rms(2,:) = rms(2,:)+ (jauto(2,i,:)-jautoavg(2,:))**2
        rms1(2,:) = rms1(2,:) +(j1auto(2,i,:)-j1autoavg(2,:))**2

  enddo

  rms(:,:) = SQRT(rms(:,:)/10)
  rms1(:,:) = SQRT(rms1(:,:)/10)


!!!!!!!!!!!calculate tau and tau2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   jjint = 0.
  jjint1 = 0.
   do isnap=0,nsnaps
         jjint(:) = jjint(:)+jautoavg(:,isnap)-jautoavg(:,(nsnaps/2))
         jjint1(:) = jjint1(:)+j1autoavg(:,isnap)
   enddo
         jjint = jjint * snaptime
         jjint1 = jjint1 * snaptime
         print *, 'average of tau of delta_J = ', jjint

!write tau and tau2 into file
  open(file='x_tau_tau2.dat',status='unknown',unit=10)
  write(10,9219)jjint(1),jjint1(1)
9219 format(1x,2g24.16)

 open(file='y_tau_tau2.dat',status='unknown',unit=10)
  write(10,9220)jjint(2),jjint1(2)
9220 format(1x,2g24.16)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


! write a file with average of autocorrelation of delta_j for mode x
   open(file='x_jautoavg.dat',status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9215)isnap*snaptime,jautoavg(1,isnap),rms(1,isnap)
  9215 format(1x,2g24.16,2g24.16)
  enddo
  close(10)

! write a file with average of autocorrelation of j_dot for mode x
  open(file='x_j1autoavg.dat',status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9216)isnap*snaptime,j1autoavg(1,isnap),rms1(1,isnap)
9216 format(1x,2g24.16,2g24.16)
  enddo
  close(10)
! write a file with average of autocorrelation of delta_j for mode y  
   open(file='y_jautoavg.dat',status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9217)isnap*snaptime,jautoavg(2,isnap),rms(2,isnap)
9217 format(1x,2g24.16,2g24.16)
  enddo
  close(10)

! write a file with average of autocorrelation of j_dot for mode y
  open(file='y_j1autoavg.dat',status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9218)isnap*snaptime,j1autoavg(2,isnap),rms1(2,isnap)
9218 format(1x,2g24.16,2g24.16)
  enddo
  close(10)

 call closeemc

 deallocate(jjarray)
 deallocate(in)
 deallocate(out)
 deallocate(jauto)
 deallocate(jjarray1)
 deallocate(in1)
 deallocate(out1)
 deallocate(j1auto)
 deallocate(jautoavg)
 deallocate(j1autoavg)
 deallocate(in2)
 deallocate(out2)
 deallocate(in3)
 deallocate(out3)
 deallocate(rms)
 deallocate(rms1)

END SUBROUTINE mc_md


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MC&MD without FFT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine mc(qs,ps)
  
  use prec
  use global_data
  use md_vv
  use emc
  use global_data
  IMPLICIT NONE

    real(r8), dimension(:,:), allocatable ::jautoavg, j1autoavg, rms,rms1
    real(r8),dimension(:,:,:), allocatable :: jauto, j1auto
    real(r8), dimension(:,:), allocatable :: jjarray, aarray, aaarray, aarray1,aaarray1
    real(r8):: num,  den, taugk, den1, chi
    real(r8), dimension(2) :: jtot, j0, j02, j1, j1tot, javg, action,laction, action0, laction0, jjint, jjint1
    real(r8), dimension(2), intent(in) :: qs, ps
    real(r8), dimension(2) :: us, vs
    real(r8), dimension(2) :: uscopy, vscopy
    integer(i12) :: isample, msample         ! sample index
    integer :: isnap, i
    real(r8) :: emccurrent   ! current value of energy as calculated by emc routines
  
  ! for output file
  character (len=20),dimension(2,20) :: correlation  ! file for correlation fn for delta_J and J_dot mode x and y

correlation(1,1) = "x_correl1.dat"  ! file for autocorrelation fn for delta_j of mode x on trial 1
correlation(1,2) = "x_correl2.dat"
correlation(1,3) = "x_correl3.dat"
correlation(1,4) = "x_correl4.dat"
correlation(1,5) = "x_correl5.dat"
correlation(1,6) = "x_correl6.dat"
correlation(1,7) = "x_correl7.dat"
correlation(1,8) = "x_correl8.dat"
correlation(1,9) = "x_correl9.dat"
correlation(1,10) = "x_correl10.dat"
correlation(1,11) = "x_dcorrel1.dat" ! file for autocorrelation fn for j_dot of mode x on trial 1
correlation(1,12) = "x_dcorrel2.dat"
correlation(1,13) = "x_dcorrel3.dat"
correlation(1,14) = "x_dcorrel4.dat"
correlation(1,15) = "x_dcorrel5.dat"
correlation(1,16) = "x_dcorrel6.dat"
correlation(1,17) = "x_dcorrel7.dat"
correlation(1,18) = "x_dcorrel8.dat"
correlation(1,19) = "x_dcorrel9.dat"
correlation(1,20) = "x_dcorrel10.dat"
correlation(2,1) = "y_correl1.dat"   ! file for autocorrelation fn for delta_j of mode y on trial 1
correlation(2,2) = "y_correl2.dat"
correlation(2,3) = "y_correl3.dat"
correlation(2,4) = "y_correl4.dat"
correlation(2,5) = "y_correl5.dat"
correlation(2,6) = "y_correl6.dat"
correlation(2,7) = "y_correl7.dat"
correlation(2,8) = "y_correl8.dat"
correlation(2,9) = "y_correl9.dat"
correlation(2,10) = "y_correl10.dat"
correlation(2,11) = "y_dcorrel1.dat"  ! file for autocorrelation fn for j_dot of mode y on trial 1
correlation(2,12) = "y_dcorrel2.dat"
correlation(2,13) = "y_dcorrel3.dat"
correlation(2,14) = "y_dcorrel4.dat"
correlation(2,15) = "y_dcorrel5.dat"
correlation(2,16) = "y_dcorrel6.dat"
correlation(2,17) = "y_dcorrel7.dat"
correlation(2,18) = "y_dcorrel8.dat"
correlation(2,19) = "y_dcorrel9.dat"
correlation(2,20) = "y_dcorrel10.dat"


  allocate( jjarray(2,0:nsnaps) )
  allocate( aarray(2,0:nsnaps) )
  allocate( aaarray(2,0:nsnaps) )
  allocate( aarray1(2,0:nsnaps) )
  allocate( aaarray1(2,0:nsnaps) )
  allocate(jautoavg(2,0:nsnaps))
  allocate(j1autoavg(2,0:nsnaps))
  allocate(rms(2,0:nsnaps))
  allocate(rms1(2,0:nsnaps))
  allocate(jauto(2,1:10,0:nsnaps))
  allocate(j1auto(2,1:10,0:nsnaps))


  us = qs
  vs = ps
  jautoavg = 0.
  j1autoavg = 0.
  rms = 0.
  rms1 = 0.

  do i=1,10

  jjarray = 0.
  j02 = 0.

  aarray = 0.
  aaarray = 0.
  aarray1 = 0.
  aaarray1 = 0.

  call setupemc( du, dv, us, vs, m, seed+ 2*i, lambda, emccurrent )
  call setupmdvv( lambda, dt )

!+++++++ equilibrate first
!
  call manystepsemc( us, vs, m, nmcequil )
!
!

!+++++++ now start sampling
  
  msample = max(nsamples/100,1)

  do isample = 1,nsamples

     if ( mod(isample-1,msample).eq.0 ) print *,' isample = ',isample

     call manystepsemc( us, vs, m, nmcsteps )

        j0(1) = (vs(1)**2 + us(1)**2*omega1**2)/(2*omega1)
        j0(2) = (vs(2)**2 + m**2*us(2)**2*omega2**2)/(2*omega2*m)
        j1(1) = (vs(1)*us(1)*(1-omega1**2+lambda*us(2)**2))/omega1
        j1(2) = (vs(2)*us(2)*(1-m*omega2**2+lambda*us(1)**2))/(omega2*m)
	print *,'us(2)=', us(2)

! copy the initial state so we can also do a trajectory with p ->  -p
     uscopy = us
     vscopy = vs
     
! normal modes at beginning of MD
!     print *,'zeroth term beginning at: ', aarray(2,0)
     action = j0
     aarray(:,0) = aarray(:,0) + action(:)
     aaarray(:,0) = aaarray(:,0) + action(:)**2
     laction = j1
     aarray1(:,0) = aarray1(:,0) +laction(:)
     aaarray1(:,0) = aaarray1(:,0) + laction(:)**2
     action0 = action
     laction0 = laction

     j02 = j02 + j0**2
     jjarray(:,0) = jjarray(:,0) + j0(:)**2

!     print *,'zeroth term: ',aarray(2,0) 
     
     do isnap=1,nsnaps
        call mdvv( us, vs, m, nmdstepsbetweensnaps )
!        print *, 'x = ', us

        jtot(1) = (vs(1)**2 + us(1)**2*omega1**2)/(2*omega1)
        jtot(2) = (vs(2)**2 + m**2*us(2)**2*omega2**2)/(2*omega2*m)
        j1tot(1) = (vs(1)*us(1)*(1-omega1**2+lambda*us(2)**2))/omega1
        j1tot(2) = (vs(2)*us(2)*(1-m*omega2**2+lambda*us(1)**2))/(omega2*m)

        jjarray(:,isnap) = jjarray(:,isnap) + j0(:)*jtot(:)
        action = jtot
        laction = j1tot
        aarray(:,isnap) = aarray(:,isnap) + action(:)
        aaarray(:,isnap) = aaarray(:,isnap) + action(:)*action0(:)
        aarray1(:,isnap) = aarray1(:,isnap) + laction(:)
        aaarray1(:,isnap) = aaarray1(:,isnap) + laction(:)*laction0(:) 
!	print *, aarray(2,isnap)
     enddo

     us = uscopy
     vs = -vscopy

! action at beginning of MD
        action(1) = (vs(1)**2 + us(1)**2*omega1**2)/(2*omega1)
        action(2) = (vs(2)**2 + m**2*us(2)**2*omega2**2)/(2*omega2*m)

     aarray(:,0) = aarray(:,0) + action(:)
     aaarray(:,0) = aaarray(:,0) + action(:)**2
     action0 = action

        laction(1) = (vs(1)*us(1)*(1-omega1**2+lambda*us(2)**2))/omega1
        laction(2) = (vs(2)*us(2)*(1-m*omega2**2+lambda*us(1)**2))/(omega2*m)

     aarray1(:,0) = aarray1(:,0) + laction(:)
     aaarray1(:,0) = aaarray1(:,0) + laction(:)**2
     laction0 = laction

     j02 = j02 + j0**2
     jjarray(:,0) = jjarray(:,0) + j0(:)**2

!     print *,'zeroth term: ',aarray(2,0)

     do isnap=1,nsnaps
        call mdvv( us, vs, m, nmdstepsbetweensnaps )

        jtot(1) = (vs(1)**2 + us(1)**2*omega1**2)/(2*omega1)
        jtot(2) = (vs(2)**2 + m**2*us(2)**2*omega2**2)/(2*omega2*m)
        j1tot(1) = (vs(1)*us(1)*(1-omega1**2+lambda*us(2)**2))/omega1
        j1tot(2) = (vs(2)*us(2)*(1-m*omega2**2+lambda*us(1)**2))/(omega2*m)

        jjarray(:,isnap) = jjarray(:,isnap) + j0(:)*jtot(:)
        action = jtot 
        laction = j1tot
        aarray(:,isnap) = aarray(:,isnap) + action(:)
        aaarray(:,isnap) = aaarray(:,isnap) + action(:)*action0(:)
        aarray1(:,isnap) = aarray1(:,isnap) + laction(:)
        aaarray1(:,isnap) = aaarray1(:,isnap) + laction(:)*laction0(:)
!	print *, aarray(2,isnap)
     enddo
     
  enddo
  jjarray = jjarray/(2*nsamples)
  j02 = j02/(2*nsamples)
  aarray = aarray/(2*nsamples)
  aaarray = aaarray/(2*nsamples)
  aarray1 = aarray1/(2*nsamples)
  aaarray1 = aaarray1/(2*nsamples)
  
 ! do isnap=0,nsnaps
!     print *, 'zeroth term: ',aaarray(2,0)
 ! enddo

!  print *,' jjarray = '
!  t0 = 0.
!  print 9222,t0,j0**2
!9222 format(1x,2g20.10)
!  do isnap=1,nsnaps
!     t = isnap*snaptime
!     print 9222,t,jjarray(isnap)
!  enddo

! write a file with autocorrelation fn for delta_J of mode x
  open(file=correlation(1,i),status='unknown',unit=10)
  den = aaarray(1,0) - aarray(1,0)**2
  do isnap=0,nsnaps
     num = aaarray(1,isnap)-aarray(1,0)*aarray(1,isnap)
     chi = num/den
     jautoavg(1,isnap) = jautoavg(1,isnap) + chi
     jauto(1,i,isnap) = chi
     write(10,9211)isnap*snaptime,chi
9211 format(1x,2g24.16)
  enddo
  close(10)

! write a file with autocorrelation fn for J_dot of mode x
  open(file=correlation(1,10+i),status='unknown',unit=10)
  den1 = aaarray1(1,0)
  do isnap=0,nsnaps
     num = aaarray1(1,isnap)
     chi = num/den1
     j1autoavg(1,isnap) = j1autoavg(1,isnap) + chi
     j1auto(1,i,isnap) = chi
     write(10,9212)isnap*snaptime,chi
9212 format(1x,2g24.16)
  enddo
  close(10)
  
  ! write a file with autocorrelation fn for delta_J of mode y
  open(file=correlation(2,i),status='unknown',unit=10)
  den = aaarray(2,0) - aarray(2,0)**2
  do isnap=0,nsnaps
     num = aaarray(2,isnap)-aarray(2,0)*aarray(2,isnap)
     chi = num/den
     jautoavg(2,isnap) = jautoavg(2,isnap) + chi
     jauto(2,i,isnap) = chi
     write(10,9213)isnap*snaptime,chi
9213 format(1x,2g24.16)
  enddo
  close(10)

! write a file with autocorrelation fn for J_dot of mode y
  open(file=correlation(2,10+i),status='unknown',unit=10)
  den1 = aaarray1(2,0)
  do isnap=0,nsnaps
     num = aaarray1(2,isnap)
     chi = num/den1
     j1autoavg(2,isnap) = j1autoavg(2,isnap) + chi
     j1auto(2,i,isnap) = chi
     write(10,9214)isnap*snaptime,chi
9214 format(1x,2g24.16)
  enddo
  close(10)

  enddo

    jautoavg = jautoavg /10
  j1autoavg =j1autoavg /10

do i=1,10
	rms(1,:) = rms(1,:) + (jauto(1,i,:)-jautoavg(1,:))**2
	rms(2,:) = rms(2,:) + (jauto(2,i,:)-jautoavg(2,:))**2
	rms1(1,:) = rms1(1,:) + (j1auto(1,i,:)-j1autoavg(1,:))**2 
	rms1(2,:) = rms1(2,:) + (j1auto(2,i,:)-j1autoavg(2,:))**2
enddo

  rms = SQRT(rms/10)
  rms1 = SQRT(rms1/10)

!!!!!!!!!!!calculate tau and tau2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   jjint = 0.
  jjint1 = 0.
   do isnap=0,nsnaps
         jjint(:) = jjint(:)+jautoavg(:,isnap)-jautoavg(:,nsnaps)
         jjint1(:) = jjint1(:)+j1autoavg(:,isnap)
   enddo
         jjint = jjint * snaptime
         jjint1 = jjint1 * snaptime
         print *, 'average of tau of delta_J = ', jjint
  open(file='x_tau_tau2.dat',status='unknown',unit=10)
  write(10,9219)jjint(1),jjint1(1)
9219 format(1x,2g24.16)
!write tau and tau2 into file
 open(file='y_tau_tau2.dat',status='unknown',unit=10)
  write(10,9220)jjint(2),jjint1(2)
9220 format(1x,2g24.16)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 
! write a file with average of autocorrelation of delta_j for mode x
   open(file='x_jautoavg.dat',status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9215)isnap*snaptime,jautoavg(1,isnap),rms(1,isnap)
  9215 format(1x,2g24.16,2g24.16)
  enddo
  close(10)

! write a file with average of autocorrelation of j_dot for mode x
  open(file='x_j1autoavg.dat',status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9216)isnap*snaptime,j1autoavg(1,isnap),rms1(1,isnap)
9216 format(1x,2g24.16,2g24.16)
  enddo
  close(10)

! write a file with average of autocorrelation of delta_j for mode y
open(file='y_jautoavg.dat',status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9217)isnap*snaptime,jautoavg(2,isnap),rms(2,isnap)
  9217 format(1x,2g24.16,2g24.16)
  enddo
  close(10)

! write a file with average of autocorrelation of j_dot for mode y
  open(file='y_j1autoavg.dat',status='unknown',unit=10)
  do isnap=0,nsnaps
     write(10,9218)isnap*snaptime,j1autoavg(2,isnap),rms1(2,isnap)
9218 format(1x,2g24.16,2g24.16)
  enddo
  close(10)

  
  call closeemc

  deallocate(jjarray)
  deallocate(aarray)
  deallocate(aaarray)
  deallocate(aarray1)
  deallocate(aaarray1)
  deallocate(jautoavg)
  deallocate(j1autoavg)
  deallocate(jauto)
  deallocate(j1auto)
  deallocate(rms)
  deallocate(rms1)


end subroutine mc
